QK Box 12.2

Chadha et al. (2019) set up a completely randomized experiment studying the effects of different soil water holding capacities (WHCs) on various characteristics of the agronomic weed Lactuca serriola. Plants were grown individually in pots allocated to one of four WHCs (100%, 75%, 50% and 25%). There were seven plants for each WHC, although we just use a subset of the data from the 100% treatment group so there was no between-subjects factor. The number of leaves on each of the seven plants in 100% WHC soil was recorded weekly for nine weeks (counts at the start of the experiment, week 0, were omitted). Time was the within-subjects (repeated measures) fixed factor and individual plants were the random subjects. With nine weeks and reasonably linear trends through time for each plant (Figure 12.2), it made sense to treat time as a continuous covariate for analysis. Although week 0 was not included in the analysis, we did not centre time for analysis, so intercepts represent the number of leaves for week 0.

Prickly lettuce or Milk thistle. Mick Keough, CC SA-BY 4.0

**Link to paper **

Chadha, A., Florentine, S., Chauhan, B. S., Long, B. & Jayasundera, M. (2019). Influence of soil moisture regimes on growth, photosynthetic capacity, leaf biochemistry and reproductive capabilities of the invasive agronomic weed; Lactuca serriola. PLoS One, 14, e0218191.

Preliminaries

First, load the required packages (afex, car, lattice, lme4, lmerTest, nlme, VCA, ez, emmeans, Rmisc, MuMIn)

Import chadha data file

chadha <- read.csv("../data/chadha.csv")
chadha

set contrasts from afex

set_sum_contrasts()

setting contr.sum globally: options(contrasts=c('contr.sum', 'contr.poly'))

select 100%WHC

chadha1 <- subset(chadha, treat=="100WHC")

make individual a factor

chadha1$plant <- factor(chadha1$plant)

plot slopes - very consistent

xyplot(noleaves~week|plant, type=c("p","r"), auto.key=T, chadha1)

Fit OLS “ancova” model

chadha1.aov <- aov(noleaves~week*plant, data=chadha1)
plot(chadha1.aov)

summary(chadha1.aov)

            Df Sum Sq Mean Sq F value  Pr(>F)    
week         1   4999    4999  310.55 < 2e-16 ***
plant        6    710     118    7.35 1.2e-05 ***
week:plant   6    110      18    1.14    0.35    
Residuals   49    789      16                    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Get Type III SS)

Anova(lm(chadha1.aov), type='III')

Anova Table (Type III tests)

Response: noleaves
            Sum Sq Df F value Pr(>F)    
(Intercept)   3568  1  221.63 <2e-16 ***
week          4999  1  310.55 <2e-16 ***
plant          131  6    1.35   0.25    
week:plant     110  6    1.14   0.35    
Residuals      789 49                   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Get correct F-ratio and P value for week (tested against week by plant)

4999.050/18.400

[1] 271.69

1-pf(4999.050/18.400, 1, 6, lower.tail = TRUE, log.p = FALSE)

[1] 3.1782e-06

no GG and HF adjustments as week is continuous with 1 df

Get variance components

chadha1.vca <- anovaMM(noleaves~week+(plant)+(week*plant), chadha1)
chadha1.vca

ANOVA-Type Estimation of Mixed Model:
--------------------------------------

    [Fixed Effects]

   int   week 
16.401  3.450 


    [Variance Components]

  Name       DF        SS         MS         VC        %Total    SD       CV[%]    
1 total      15.101789                       26.537662 100       5.151472 15.308619
2 plant      6         710.095238 118.349206 10.401893 39.196719 3.225197 9.584311 
3 week:plant 6         110.4      18.4       0.038377  0.144613  0.1959   0.582156 
4 error      49        788.772222 16.097392  16.097392 60.658668 4.012156 11.922915

Mean: 33.651 (N = 63) 

Experimental Design: balanced  |  Method: ANOVA

VCAinference(chadha1.vca, alpha=0.05, VarVC=TRUE, excludeNeg=FALSE, constrainCI=FALSE)

Inference from Mixed Model Fit
------------------------------

> VCA Result:
-------------

    [Fixed Effects]

   int   week 
16.401  3.450 


    [Variance Components]

  Name       DF      SS       MS       VC      %Total  SD     CV[%]   Var(VC)
1 total      15.1018                   26.5377 100     5.1515 15.3086        
2 plant      6       710.0952 118.3492 10.4019 39.1967 3.2252 9.5843  77.6065
3 week:plant 6       110.4    18.4     0.0384  0.1446  0.1959 0.5822  0.0343 
4 error      49      788.7722 16.0974  16.0974 60.6587 4.0122 11.9229 10.5766

Mean: 33.651 (N = 63) 

Experimental Design: balanced  |  Method: ANOVA


> VC:
-----
           Estimate  CI LCL  CI UCL One-Sided LCL One-Sided UCL
total       26.5377 14.5063 63.3391       15.9486       54.6613
plant       10.4019 -6.8643 27.6681       -4.0884       24.8922
week:plant   0.0384 -0.3245  0.4013       -0.2662        0.3429
error       16.0974 11.2325 24.9968       11.8901       23.2468

> SD:
-----
           Estimate  CI LCL CI UCL One-Sided LCL One-Sided UCL
total        5.1515  3.8087 7.9586        3.9936        7.3933
plant        3.2252 -2.6200 5.2600       -2.0220        4.9892
week:plant   0.1959 -0.5697 0.6335       -0.5159        0.5856
error        4.0122  3.3515 4.9997        3.4482        4.8215

> CV[%]:
--------
           Estimate  CI LCL  CI UCL One-Sided LCL One-Sided UCL
total       15.3086 11.3183 23.6505       11.8677       21.9707
plant        9.5843 -7.7858 15.6313       -6.0087       14.8264
week:plant   0.5822 -1.6929  1.8825       -1.5332        1.7403
error       11.9229  9.9596 14.8575       10.2470       14.3280


95% Confidence Level  
SAS PROC MIXED method used for computing CIs 

Drop interaction

chadha2.aov <- aov(noleaves~week+plant, data=chadha1)
Anova(lm(chadha2.aov), type='III')

Anova Table (Type III tests)

Response: noleaves
            Sum Sq Df F value Pr(>F)    
(Intercept)   3568  1  218.22 <2e-16 ***
week          4999  1  305.78 <2e-16 ***
plant          710  6    7.24  1e-05 ***
Residuals      899 55                   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

get new variance components

chadha2.vca <- anovaMM(noleaves~week+(plant), chadha1)
chadha2.vca

ANOVA-Type Estimation of Mixed Model:
--------------------------------------

    [Fixed Effects]

   int   week 
16.401  3.450 


    [Variance Components]

  Name  DF        SS         MS         VC        %Total    SD       CV[%]    
1 total 23.462944                       27.681988 100       5.261368 15.635196
2 plant 6         710.095238 118.349206 11.333402 40.941432 3.366512 10.004257
3 error 55        899.172222 16.348586  16.348586 59.058568 4.043338 12.015581

Mean: 33.651 (N = 63) 

Experimental Design: balanced  |  Method: ANOVA

VCAinference(chadha2.vca, alpha=0.05, VarVC=TRUE, excludeNeg=FALSE, constrainCI=FALSE)

Inference from Mixed Model Fit
------------------------------

> VCA Result:
-------------

    [Fixed Effects]

   int   week 
16.401  3.450 


    [Variance Components]

  Name  DF      SS       MS       VC      %Total  SD     CV[%]   Var(VC)
1 total 23.4629                   27.682  100     5.2614 15.6352        
2 plant 6       710.0952 118.3492 11.3334 40.9414 3.3665 10.0043 57.76  
3 error 55      899.1722 16.3486  16.3486 59.0586 4.0433 12.0156 9.7191 

Mean: 33.651 (N = 63) 

Experimental Design: balanced  |  Method: ANOVA


> VC:
-----
      Estimate  CI LCL CI UCL One-Sided LCL One-Sided UCL
total   27.682 16.7947 54.045       18.1687        48.323
plant   11.333 -3.5623 26.229       -1.1675        23.834
error   16.349 11.6201 24.704       12.2651        23.081

> SD:
-----
      Estimate  CI LCL CI UCL One-Sided LCL One-Sided UCL
total   5.2614  4.0981 7.3516        4.2625        6.9515
plant   3.3665 -1.8874 5.1214       -1.0805        4.8820
error   4.0433  3.4088 4.9703        3.5022        4.8042

> CV[%]:
--------
      Estimate  CI LCL CI UCL One-Sided LCL One-Sided UCL
total   15.635 12.1784 21.847       12.6668        20.658
plant   10.004 -5.6088 15.219       -3.2109        14.508
error   12.016 10.1300 14.770       10.4073        14.277


95% Confidence Level  
SAS PROC MIXED method used for computing CIs 

Now mixed effects modelcomparing random slopes and random intercepts using ML

chadha.lmer1 <- lmer(noleaves~week + (week|plant), REML=FALSE, chadha1)

boundary (singular) fit: see help('isSingular')

 # singular fit so set correlation between random effects to zero
chadha.lmer1 <- lmer(noleaves~week + (week||plant), REML=FALSE, chadha1)
chadha.lmer2 <- lmer(noleaves~week + (1|plant), REML=FALSE, chadha1)

Compare models

anova(chadha.lmer1, chadha.lmer2)

Data: chadha1
Models:
chadha.lmer2: noleaves ~ week + (1 | plant)
chadha.lmer1: noleaves ~ week + (week || plant)
             npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
chadha.lmer2    4 375 383   -183      367                    
chadha.lmer1    5 376 386   -183      366  0.83  1       0.36

AICc(chadha.lmer1, chadha.lmer2)

Focus on random intercept model - refit with REML

chadha.lmer3 <- lmer(noleaves~week + (1|plant), REML=TRUE, chadha1)
summary(chadha.lmer3, ddf="Kenward-Roger")

Linear mixed model fit by REML. t-tests use Kenward-Roger's method ['lmerModLmerTest']
Formula: noleaves ~ week + (1 | plant)
   Data: chadha1

REML criterion at convergence: 365.6

Scaled residuals: 
   Min     1Q Median     3Q    Max 
-1.923 -0.509  0.086  0.320  3.927 

Random effects:
 Groups   Name        Variance Std.Dev.
 plant    (Intercept) 11.3     3.37    
 Residual             16.3     4.04    
Number of obs: 63, groups:  plant, 7

Fixed effects:
            Estimate Std. Error     df t value Pr(>|t|)    
(Intercept)   16.401      1.689 13.433    9.71  1.9e-07 ***
week           3.450      0.197 55.000   17.49  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
     (Intr)
week -0.584

anova(chadha.lmer3, type=3, ddf="Kenward-Roger")

Type III Analysis of Variance Table with Kenward-Roger's method
     Sum Sq Mean Sq NumDF DenDF F value Pr(>F)    
week   4999    4999     1    55     306 <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Get variance components

chadha.ci3 <- confint.merMod(chadha.lmer3, oldNames=FALSE)

Computing profile confidence intervals ...

chadha.vc3 <- (chadha.ci3)^2
print(chadha.vc3)

                        2.5 %  97.5 %
sd_(Intercept)|plant   2.7657  38.240
sigma                 11.3276  23.828
(Intercept)          170.1183 390.404
week                   9.3643  14.745

Random slopes output for comparison - refit using REML

chadha.lmer4 <- lmer(noleaves~week + (week||plant), REML=TRUE, chadha1)
summary(chadha.lmer4, ddf="Kenward-Roger")

Linear mixed model fit by REML. t-tests use Kenward-Roger's method ['lmerModLmerTest']
Formula: noleaves ~ week + (week || plant)
   Data: chadha1

REML criterion at convergence: 364.6

Scaled residuals: 
   Min     1Q Median     3Q    Max 
-2.017 -0.521  0.077  0.375  3.925 

Random effects:
 Groups   Name        Variance Std.Dev.
 plant    (Intercept)  6.797   2.607   
 plant.1  week         0.141   0.375   
 Residual             15.644   3.955   
Number of obs: 63, groups:  plant, 7

Fixed effects:
            Estimate Std. Error    df t value Pr(>|t|)    
(Intercept)    16.40       1.47 10.28    11.2  4.4e-07 ***
week            3.45       0.24 10.28    14.4  3.8e-08 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
     (Intr)
week -0.530

anova(chadha.lmer4, type=3, ddf="Kenward-Roger")

Type III Analysis of Variance Table with Kenward-Roger's method
     Sum Sq Mean Sq NumDF DenDF F value  Pr(>F)    
week   3245    3245     1  10.3     207 3.8e-08 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Check whether AR(1) covariance structure would improve the fit using nlme

chadha.lme1 <- lme(noleaves~week, random=~1|plant, method="ML", chadha1)
chadha.lme2 <- lme(noleaves~week, random=~1|plant, method="ML", correlation=corAR1(form=~1|plant), chadha1)
anova(chadha.lme1, chadha.lme2)

AICc(chadha.lme1, chadha.lme2)

no improvement with AR(1)

chadha <- lm(noleaves~week, chadha1)
summary(chadha)

Call:
lm(formula = noleaves ~ week, data = chadha1)

Residuals:
   Min     1Q Median     3Q    Max 
-9.101 -3.151  0.249  2.949 14.549 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   16.401      1.410    11.6   <2e-16 ***
week           3.450      0.251    13.8   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 5.14 on 61 degrees of freedom
Multiple R-squared:  0.756, Adjusted R-squared:  0.752 
F-statistic:  189 on 1 and 61 DF,  p-value: <2e-16